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1 Introduction 



Models of population dynamics have been the subject of much study in bi- 



ology. The essay of Malthus (1798) first suggested exponential growth. The 



idea of quantitative hmits on growth is usually attributed to Verhulst (1838), 
who pioneered the use of the logistic equation (LE), and to Lotka (1920) and 



Volterra (1926), who considered models of predators and prey. Typically these 



early studies of population dynamics regards the units of interest as the num- 
ber of individuals in populations of one or more species, without regard to 
genetic diversity within a population. 

Given an initial focus on large population sizes, this early work focused on de- 
terministic equations with continuous rather than discrete variables. A related 
approach is the logistic map (LM) which is also deterministic and continuous. 



but treats time as a series of discrete generations (May, 1976). Continuous 



population equations like this can be extended to include external environ- 
mental fluctuations. For continuous time, this results in a continuous stochas- 
tic logistic equation (SLE) used to model populations: typically for ecological 
studies. While it is possible to extend continuous state models to incorporate 
genetics, they are limited by the fact that reproduction and mutation have 
an intrinsically random nature, requiring a statistical treatment of a discrete 
jump process. 



The discrete logistic equation (DLE) was introduced by FeUer (1939). This 
is a continuous time Markov process with discrete states. It describes the 
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probabilistic time-evolution of a discrete rather than continuous population, in 
which the rates of birth or death are linear or quadratic in the total population 
number. This model and variations of it have been the subject of study by a 
number of researchers, mostly in the context of calculating extinction times 



(Nasell 


2001 


Matis et al. 


2003 



Doering et al. , 2005). One of the principle virtues of the model is that it 
can treat demographic noise: the inevitable stochasticity intrinsic to discrete 
random events such as birth and death. 



A particular special case of stochastic logistic growth is obtained when density 
dependence occurs only through asymptotically decreasing birth rates, such 
that a hard upper limit in population size (known as the maximum carry- 



ing capacity) can never be exceeded (Dushoff, 2000). This model originally 



arose in mathematical epidemiology as the closed stochastic SIS model (Weiss 



and Dishon, 1971; Nasell, 1996); a simple model of endemic disease in which 
the maximum carrying capacity is the (constant) size of the host population. 
It is also possible to construct a discrete logistic equation that has density 
dependent death rates, which is the approach we use here. 



Any continuous-time discrete-state Markov process of this type can be canoni- 
cally described by a birth-death master equation and such equations are widely 



used in chemistry and physics (Gardiner, 2004). An approximate way to treat 
Master equations - introduced by Kramers and Moyal - is to replace the dis- 
crete variables (like population size) by continuous variables undergoing a dif- 
fusion process( Kramers, 1940; MoyaH[l949 ). With truncation, this results in a 
Fokker-Planck equation. After transforming to a continuous variable stochas- 



tic equation, this is simply equivalent to the SLE (Feller, 1951). 



This approximation results in a form of the SLE which can describe the ef- 
fects of demographic noise in the limit of large population numbers. Thus, the 
Kramers-Moyal approximation to the discrete logistic equation is similar to 
the stochastic logistic equation for environmental noise, but with a different 
interpretation of the noise source. It has the drawback that it is not accurate 



for small populations. This leads to exponentially large errors(Gaveau et al. 



1996) for important problems like extinction times, which necessarily involve 
a small population. Nevertheless, many pubhshed works in mathematical ecol- 
ogy and epidemiology that address demographic stochasticity have used the 
Kramers-Moyal or even more serious types of approximation. 



At the same time, quantitative statistical models of the genetic evolution 
of populations are based on the pioneering work of population geneticists 



Fisher (Fisher, 1918, 1930), Wright (1930) and Haldane (1932), whose com 



bined work has laid the foundation for subsequent research in this area. How- 
ever, many important questions are inaccessible by current techniques, be- 
cause most population genetic theory is based on the ideahzed Wright-Fisher 
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model of population and the coalescent (Kingman, 1982), neither of which 



handle small randomly fluctuating populations. Thus, the flelds of mathemat- 
ical ecology and epidemiology, where questions about extinction, population 



viability (Boyce 1992 Lande 1993 Foley, 1994; Engen and Saether 2000 



Hakoyama and Iwasa, 2000; Drake, 2006; Cairns et al., 2007) and short-term 



population dynamics are important, have developed quite independently from 
theoretical population genetics. 



Recently there has been renewed interest in unifying these flelds from two 
different directions. Firstly researchers interested in the study of viral evolu- 
tionary dynamics have long realized that epidemiological dynamics and pop- 



ulation genetics have overlapping time scales in viruses (Pybus et al., 2001 



2003; Kelly et al., 2003; GrenfeU et al. , 2004). Secondly and more recently 



there have been a number of attempts to develop general theoretical results 
aimed at extending classical population genetics concepts and results to the 
analysis of self-regulating stochastic population models. Notable examples in- 
clude calculation of the effective population size of fluctuating populations 



that have both environmental and demographic stochasticity (Engen et al 



2007); development of a branching-process analogue of the stochastic logistic 



model with linear birth rates and density-dependent death rates (Lambert 



2005) and its apphcation to the analysis of flxation probabilities in statisti- 



cally varying populations (Lambert, 2006; Champagnat and Lambert, 2007); 



and analysis of the one-locus two-aUele flxation probability in closed stochastic 



SIS model (Parsons and Quince, 2007). 



Similar underlying mathematical issues to those faced by biologists have been 
known in the physical sciences for some time. Efficient techniques for solving 



master equations directly by simulation were introduced by GiUespie (1977) 



An exact transformation of master equations to a continuous Fokker-Planck 



equation and hence a stochastic form was introduced by Gardiner (2004) 



and is known as the Poisson representation. This approach differs from the 
commonly used Kramers-Moyal approximation in that it is exact for small 
numbers, provided boundary terms vanish, and hence can be used to rehably 
calculate extinction times (among other things). The more recently developed 



stochastic-gauge Poisson representation(Deuar and Drummond, 2002; Drum 



mond and Deuar, 2003; Drummond, 2004) eliminates boundary term errors 



that can occur in the Poisson method, thus making the technique more gen- 
erally useful. 



In this paper, we apply these theoretical techniques to the problem of calcu- 
lating and simulating extinction times in the discrete logistic equation with 
environmental noise added. We consider a logistic model that includes demo- 
graphic noise due to birth, death, density-dependent intraspeciflc competition 
and a simple model of environmental noise. We obtain exact analytic results in 
this equation which unifles the DLE and SLE models. For pure demographic 
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noise, we have verified that our exact results agree precisely with a Monte- 
Carlo simulation of the master equation, showing that the discrete jump and 
continuous noise approaches are exactly equivalent with these techniques. We 
also compare our results with direct simulations of stochastic equations and 
with steepest-descent approximations. These provide useful analytic expres- 
sions, valid in the limit of large populations. 

We present asymptotic results demonstrating that while large carrying capac- 
ity increases extinction times exponentially, it is not the only factor involved. 
There is also an exponential dependence on the reproductive ratio. For a given 
growth rate and carrying capacity, populations are exponentially shorter-lived 
at small reproductive ratio, as i? ^ 1 . 

In a following paper, we will extend these results to show that spatial connect- 
edness provides a robustness against global extinction of a population. We do 
this by describing the relationship between the rate of extinction and the size 
and connectedness of the metapopulation. This is analysed quantitatively us- 
ing a combination of analytic and numerical techniques. Extensions to include 
genetics will be treated elsewhere. 



2 Methods 



In this paper we will revisit and extend the logistic model of population 
dynamics. In its original form this model simply assumed that while small 
populations may grow exponentially, this growth saturates due to nonlinear 
effects such as density-dependent competition for resources. Thus, for a time- 
dependent population (t) of a given species, 

dN 

^=N{g-cN), (1) 



where g is the initial growth rate, c describes competition, and = g/c 
is the deterministic carrying capacity. This is known as the logistic equation 



(Verhulst, 1838). If the initial population values are small, it has the simple 



property of an initial growth followed by a gradual approach to the steady- 
state, in which Nf.. Without the competition term there is an unreahstic 
exponential growth that fails to account for the inevitable competition for 
hmited resources that occurs in all natural populations. 
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2.1 Environmental and demographic stochasticity 



Nothing is more certain than uncertainty. All populations experience random 
fluctuations in numbers. These fluctuations are intrinsic to the random timing 
of population changing events, together with the granularity of populations 
that are made up of a finite number of individuals. To describe this real- 
ity, it is necessary to use master equations that deal with discrete outcomes, 
rather than continuous rate equations. We wish to analyse this demographic 
stochasticity, as a measure of minimum, irreducible randomness which cannot 
be ehminated from real population dynamics. In nature, there is also ran- 
domness in the birth and death rates, caused by the chaotic fluctuations in 
any real environment. Obviously, this environmental stochasticity represents 
a large part of observed stochasticity, and should not be ignored. 

Traditional deterministic population equations may include environmental 
noise. However, they frequently ignore demographic stochasticity, even though 
this is essential to the analysis of extinction rates, where the discrete differ- 
ence between zero and one is important. Conversely, previous treatments of the 
discrete stochastic logistic model have generally ignored external fluctuations. 
An environmental fluctuation may change the death rate in a population, but 
at close to zero population the demographic statistical fluctuations become 
more important. These determine the minimum survivable population, before 
extinction becomes a near certainty. 

The standard treatment of genetic drift also includes stochastic fluctuations 
in the frequencies of different genotypes, but ignores such fluctuations in the 
overall number of individuals in the population. Our model can be extended 
to include both effects and thus provides a platform for further development 



in this new area (Parsons and Quince, 2007) 



2.2 Master Equations 



One can show that there is a stochastic, or random equation that corresponds 
exactly - both in terms of averages and fluctuations - to the types of process we 
are considering here. Such equations were first considered in the context of the 
logistic model, by Feller. Generically, they are known as birth-death master 
equations. We first note that all of the number fiuctuations we are interested 
in can be understood by combining the genotype and location labels together 
into one abstract species X^. 

Here, the combined label z = {j, x} indicates genotype j at location x. Thus, 
there are a set of discrete position labels x, each indicating a ceU or lattice 
site within which there is relatively strong mixing. For convenience, we label 
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these abstract combinations of location and species in this section simply as 
Xi, X2 , and so on. 



To derive the stochastic equations, first consider the most general kinetic equa- 
tion we are interested in. As a simple model for stochastic birth/death events, 
we consider a binary reaction 

kit) 

EiXi + 82X2 83X3 + £4^4 , (2) 



which describes ei individuals of type Xi combining with S2 individuals of type 
X2, to give £3 individuals of type X3 and £4 individuals of type X4. The rate 
for the process is R{t), which may be time-dependent or stochastic. There 
is a birth-death master equation whose solution gives the time-evolution of 
probabilities -P(N), for finding the total number of individuals of each type 
equal to N = {Ni,..,Nd) ■ Here D is the genetic/spatial dimension, equal 
to the number of distinct genotypes multiphed by the number of lattice sites 
under consideration. 

This master equation is obtained by summing over all processes, where each 
binary reaction generates a term of the form: 

(N) = -k{t)Nl^N^'P (N) 

+ k{t){N[Y' {N'2f' P{^') (3) 



Here, N' = (A^i -|- Si, iV2 + £2, — £3, -^4 — £4, ••, N^), with obvious modifica- 
tions in cases of identical labels. This is difficult to solve as it stands, due 
to the large dimensionality of the space of distributions over integer popu- 
lations. This complexity increases exponentially with the dimensionality 
which means that both the vector P (N) and the matrix that describe its 
microscopic changes, rapidly become too large to be stored or numerically 
analysed in the memory of any computer. 

However, there are several techniques for solving these equations which make 
use of mappings into probabilistic Markov processes or stochastic equations. 
These do not have exponential complexity, but rather trade off extremely large 
matrices in favour of sampling over a probabilistic distribution. While this 
introduces sampling errors, the amount of error can be controlled and reduced 
by increasing the number of samples. In some cases, exact or approximate 
analytic solutions can be obtained without requiring simulations. 

Two methods we will use here are the direct method of Gillespie and the 
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Poisson representation. 



2.3 Direct Gillespie Simulation 



The direct method of simulating stochastic rate equations was originally ap- 
phed to the problem of simulating chemical kinetics by Gillespie (1977). Ap- 
phed to population dynamics this technique is essentially an in-silico direct 
model of the population. In this method a single event is one reaction process, 
each event adding or subtracting an individual after a random time interval 
whose mean is the inverse of the reaction rate. 



The Gillespie algorithm produces a single forward-time reahzation of the 
stochastic process of interest. In order to obtain statistics on the average time- 
evolution of the process, many replicates are simulated, each starting with a 
different random number seed so that a sample of realisations is obtained that 
represent the range of outcomes. An individual Gillespie simulation proceeds 
one event at a time. The events considered in this paper are birth, death, and 
death by density-dependent competition, although extensions to more general 
models are easily made. 

The Gillespie method is straightforward to implement and provides a use- 
ful base-line for comparisons, especially when demographic noise is the main 
source of fluctuations. Its main limitations are that it cannot be used to model 
large populations due to its computationally intensive approach, and the ad- 
dition of environmental noise is nontrivial. 



2.4 Poisson Representation 



Whfle direct simulations are extremely useful computational tool, they are 
not directly amenable to analytic solutions, nor are they immediately suitable 
for describing a unifled theory that combines demographic and environmental 
noise. An efficient alternative technique for this purpose is to use the Pois- 



son representation (Chaturvedi et al. 


1977 


Chaturvedi and Gardiner 


1978 




Gardiner, 



Gardiner and Chaturvedi, 1977 



2004) - an exact expansion over 



Poisson distributions - to map these equations into a stochastic differential 
equation for a continuous Poisson amplitude. 



This proceeds by expanding P (N) in a basis of elementary multivariate Pois- 
son distributions, so that: 



P(N) = Jd[x\f{x)P{N;x) , 



(4) 
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where P (N |x) is the Poisson distribution, with: 



P(N;f) = n 



m 



(5) 



Each elementary process generates an equivalent Fokker-Planck equation (where 

di = d/dxi): 



d 



fit, x) = k{t) [(1 - d^r (1 - d,r - (1 - (i - ^2)^^] 



(6) 



After summing over all elementary processes, one obtains the total Fokker- 
Planck equation: 



d 
di 



/(t,x) = ^9. 



-A:(t,x) + ^^9,A,(t,x) 



/ {t, x) 



(7) 



The utility of the Poisson representation is that it can be used to directly 
calculate the factorial moments and correlations, through the equivalences: 

{x""') = {N{N -1)..{N -n + 1)) . (8) 



2.5 Stochastic realizations 



The Fokker-Planck equation can also be readily converted into an Ito stochas- 
tic differential equation for direct numerical simulation, where: 



dxl = Ai (t, x) + ^ B,j (t, x) dWj (t) 



(9) 



Here, Dij = Y^kBikBjk , and the terms dWi are independent real Gaussian 
noise terms with 



{dWidWj) = 6ijdt 



(10) 



It is necessary that Dij is positive-definite. This may require introducing a 
complex Poisson variable Xj, together with additional stochastic gauge terms to 



ensure stability. Relevant technical details are presented elsewhere (Drummond 



2004). The Ito form of stochastic equation corresponds most directly to a 
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standard Fokker-Planck equation form, which is generated directly from the 
Poisson representation of demographic noise. This describes a forward step in 
which the multiphcative noise term is evaluated at the current time. 

There is another commonly used form of stochastic equation. The Stratonovich 
equation describes a time-symmetric step, which is the broad-band limit of a 
finite band-width physical noise. This is the correct form in which to introduce 
environmental noise sources originating from time-dependent rate fiuctuations 
R{t). The two forms can be converted from one to the other, since the equiv- 
alent Stratonovich equation has the form: 

dxf = A {t, x)dt + Y^ Bij {t, x) dWj (t) (11) 
j 

where 

A {t, x) = A {t, x) - ^ E Bkj {t, x) dkBij {t, x) . (12) 

j 

These become formally identical if the noise term is independent of the phase- 
space variable x. This case is termed additive noise or constant diffusion. We 
will make use of a variable change to achieve constant diffusion behaviour in 
the analytic calculations. 



2.6 Population processes on individuals 

We can categorize biologically relevant kinetic processes into the following list, 
in which we give the simplest stochastic equation in each case. 

2. 6. 1 Transformation 

A transformation is a unary reaction in which one species changes to another 
at a constant rate, either through mutation or spatial motion (recall that the 
subscripts describe both physical and genetic space). In practise, this may 
be catalysed by other organisms or molecules, or involve additional precursor 
molecules. As long as these additional species have a constant concentration, 
they can be neglected, leading to the reaction: 

^1 X2 (13) 
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• Fokker-Planck equation: 

p-{t,x) = a[d^-d,]xj{t,x) . (14) 

• Stochastic equation (noise-free): 

dxi = —axidt 

dx2 = axidt . (15) 

As an example, the death of an organism will be treated in a simplified way, 
where the final product is neglected: 

X ^ (16) 



• Fokker-Planck equation: 

j^f{t,x)^ad,xf{t,x) . (17) 

• Stochastic equation (noise-free): 

dx — —ax dt . (18) 

None of these processes involve stochastic terms, unless the rates themselves 
fluctuate. This implies that, at constant rates, transformation processes map 
one Poisson distribution into another Poisson distribution. 

2.6.2 Birth 

Consider birth: a process in which a pre-existing organism Xi generates two 
further organisms X2 and X^. A special case of this occurs when all of the 
species involved are identical. 

b 

X 2X (19) 

• Fokker-Planck equation: 

^^f(t,x)^b[-d, + dl\xf{t,x) . (20) 

• Ito stochastic equation: 

dx = bxdt + V2bxdW (t) (21) 
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The fact that there is noise in these equations simply implies that birth is 
intrinsically noisy, in the sense that it increases local fluctuations above the 
Poisson level. 



2. 6. 3 Competition 

Next, consider competition: a process in which two pre-existing organisms Xi 
and X2 compete (usually for resources), reducing the total number to just one, 
X3. As in birth, a special cases of this occurs when all of the species involved 
are identical: 



2X 



X 



(22) 



Fokker-Planck equation: 

d 



dt 



f {t, x) 



a - di 



Ito stochastic equation: 

dx = —cx'^dt + ixV2cdW{t) 



(23) 



(24) 



Competition tends to reduce fluctuations below the Poisson level. This results 
in complex Poisson amplitudes in general, unless the fluctuations of another 
process keep the overall variance real. This equation in its present form is un- 
stable, and would require additional stochastic gauges to be treated correctly 



(Drummond, 2004). This issue does not arise in the stochastic logistic model. 



due to the additional positive diffusion from the birth term. 



3 Discrete logistic model 



Consider the discrete logistic model of a birth process X — > 2X, together with 
competition 2X X, and a death process X ^ 0. This type of random 
process can be represented kinetically as: 



X 



X 2X 



(25) 
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2X ^ X 



For small populations the resulting net growth rate is g = h — a, which cor- 
responds to the Malthusian fitness of type X. However, we have a model of 



fitness which is also density dependent and stochastic (Feller, 1939). In this 
model the density-dependence come from the death by competition. It is also 
possible to have density-dependent birth-rates, which gives an extra parameter 
at the expense of more complicated stochastic behaviour. 

Apart from the growth rate g which sets the relevant time-scales, this process is 
described by two dimensionless parameters; Nc = g/cis the carrying capacity, 
and R = b/a is the birth-death ratio, also called the reproductive ratio. 



3.1 Master Equation 



The master equation corresponding to the dynamical behaviour in Eq (25) is 



^p(N) = -{a + b + cN) NP (N) 

+ b{N -1)P{N -1) . 
+ c{N + lfP{N + l) 

+ a{N + 1)P{N + 1) (26) 



This is a set of coupled differential equations which can be treated directly 
in a single-species case, provided the total population is not too large. How- 
ever, such techniques are difficult to extend to multiple species problems, due 
to exponentially increased complexity. Therefore, in this paper we focus on 
techniques which can be readily scaled to treat more complex multi-species 
cases. 

As an example, we show a graph in Figure 1 of a single Gillespie-type realiza- 
tion of this jump process, leading to an extinction, together with the average 
over many realizations. Note that in this picture the population is always 
discrete, with integer jumps at random times. 
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12 14 16 18 20 



Figure 1. Direct simulation results using Gillespie algorithm with parameters 6 = 2, 
a = 1, c = 0.2 . This corresponds to Nc = 5, R = 2, g = 1. The initial condition 
of A'o = 8 was sampled from a Poissonian with mean of xq = 5. The solid line is a 
single stochastic realization, showing integer jump behaviour. The dotted line is an 
average of 10000 realizations, showing the exponential decline in average population 
leading to extinction. 



3.2 Poisson Equation 



The master equation (26) corresponds exactly to the following Fokker-Planck 



equation in the Poisson representation for / (x) : 



d 
dt 



fix) 



d_ 
dx 



-gx + cx^ + — (bx — cx^ 
ox 



In terms of the abstract notation of Eq (|9 



fix) 



(27) 



A{x) = gx — cx^ 
D{x) = 2{bx-cx^) 



(28) 



which implies there is a corresponding Ito stochastic differential equation: 

dx^ = (gx - cx^) dt + ^2 (bx - cx^)dW , (29) 
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where: 



{dW^) = dt . (30) 



This transforms the discrete master equation given above to a completely 
equivalent Ito stochastic differential equation. The equation is valid in the 
domain < a; < c/h , in which region the noise terms are all real. As the noise 
term vanishes at the boundaries, it is appropriate to use reflecting boundary 
conditions. The boundary at a; = is absorbing; once it is reached, the variable 
X remains at zero. This corresponds to extinction of the population on this 
trajectory. 



3.3 Numerical Solutions 



In Figure 2, we show a single stochastic reahzation of this equation, together 
with a time average. Note that here the stochastic variable x is continuous; it 
is the mean of one of the Poisson distributions used in the expansion. Even 
though there is no direct correspondence between each realization in Figures 1 
and 2, they have identical means. More general observables can be transformed 
between these stochastic methods using the factorial moment equivalences. 



In these computer simulations, we use an imphcit central difference method 



with the time-symmetric Stratonovich form of Eq (11), including the trans 



formed drift term of Eq (12), again supplemented by absorbing boundary 
conditions at a; = 0: 

dx^ ={[g + c]x- 6/2 - cx^) dt + ^2 (bx - cx^)dW , (31) 



Points to notice about this approach are: 

• The stochastic equation is for the mean of an 'equivalent' Poisson distribu- 
tion 

• P (N) must be reconstructed from an average over these Poisson distribu- 
tions 

• Even if g = there is noise if 6 7^ 0; births cannot cancel deaths exactly. 

• This intrinsic random noise occurs in addition to any environmental fluctu- 
ations 

• The relative noise size scales as of the deterministic terms in the 
equation 

• There is no noise if a; = 0; this means extinction, as shown by the Ito form. 
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Figure 2. Poisson simulation results using Stratonovich equations with parameters 
6 = 2, a = 1, c = 0.2 or i? = 2, A^c = 5 as in Figure 1. The initial condition was 
a Poissonian with mean of xq = 5. The solid line is a single stochastic realization, 
showing continuous stochastic behaviour. The dotted line is an average of 10000 
realizations, showing the exponential decline in average population leading to ex- 
tinction, just as in the Gillespie result. Integration step-size utilized was dt = 0.04, 
RK4 Runge-Kutta algorithm. Average extinction time T = 10.04zb0.27 using Eq([39|). 



3.4 Population statistics 



In order to demonstrate an important property of the model, consider the time- 
evolution of moments, without competition so that c = 0. In the stochastic 



case, moments can be calculated, using standard Ito variable change (Gar- 



diner, 2004) rules: 



{dx") = {{x + dxf) - (x") 

= n{dxx''~^) + ^ii!^^(rfxV-2) 

= ng{x'^)dt + n{n-l)b{x''-^) . (32) 



This means that, even if = 0, the population distribution does change statis- 
tically. In classical population dynamics, one would have a stable equilibrium 
with a constant population. However, demographic noise modifies this sub- 
stantially. If we start from a Poissonian with amplitude Xq, then the factorial 
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moments increase with time even though the mean is constant: 



{x^)=Xo 
{x^)=xl + 2btxo 
(x^) = xl + 6 {btf xo 
{x'')=xt + 2A{btfxo . 



(33) 



Thus, for example, the standard deviation a = y {N^ — N"^) can now be cal- 
culated to give: 



This shows that the initial Poissonian variance (equal to the population mean) 
increases hnearly with time. Some populations become large, others become 
extinct - even if g, the growth rate, is zero. 



4 Demographic extinction times 

We wish to analyse the extinction rate of a single isolated population under- 
going discrete logistic population dynamics with competition included. We 
will assume an initial Poissonian distribution with mean xq, as this allows us 
to include some of the effects of initial population fluctuations in biologically 
relevant cases. 

There are several equivalent methods to carry out this calculation. One method 
is a direct calculation from the master equation, which gives a result without 
any environmental noise effects. A direct Gillespie simulation is also feasi- 
ble. We will use this as a reference calculation, since the technique has the 
advantage that it can be readily extended to more complex cases. Because 
the simulations are relatively time-consuming for large populations and long 
extinction times, alternative approaches are desirable. 

To obtain analytic results, we will use the Poisson equation. This has both 
exact and asymptotic approximate solutions for extinction times. In the next 
section we show that the approach can be readily extended to include envi- 
ronmental noise. 




(34) 
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4-1 Exact extinction times 



Calculating first-passage times for diffusion or Fokker-Planck processes is a 
well known problem in physics, and there is an extensive literature. Here there 
is an interesting subtlety. Even when the Poisson variable is non-vanishing so 
that a; > , there is a fraction of size in the corresponding sample of 
real populations that is already extinct. To take account of this, the standard 
first-passage time calculation needs modifications. 

We start by defining the rate of extinction at time t given an initial Gaussian 
with mean Xo as 7l{t\xQ). This is given by the rate of change of the probability 
for extinction S{t\xo) given an initial Poisson distribution at xq: 



n{t\xo) = ^S{t\xo) 



dt 
d 

dt 



J f {t,x\xo)e-''dx . (35) 



The average time to extinction T given an initial Poisson distribution at xq 
involves a time-integral over the extinction rate: 



T(xo) = j tn{t\xo)dt 



= j t-£{t\xo)dt 





oo 



^- J t—A{t,\xo)dt . (36) 



Here we have introduced an alive probability defined as 

A{t,\xo)^l-£{t\xo) 

f{t,x\xo)\l-e~'']dx. (37) 







The absorbing state at x = means that the distribution decays to a delta 
function at long times: / {oo,x\xo) = S{x), so that: 

lim tA(t, \xo) = 0. 

t— »oo 
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This last result allows a further simplification from integration by parts, so 
that: 



oo 

T(xo) = J A{t, \xo) dt . 



(38) 



Since / (t,x\xo) is a normahzed, probabilistic propagator function in the Pois- 
son representation, this leads to a simple prescription for calculating extinction 
times in a numerical simulation: 



T{xo) 



(39) 



Analytically, however, we can do better, and even find an exact solution just 
requiring numerical integration. ^ is a linear functional of the Poisson distri- 
bution /(t, xlxo), and as such must satisfy the well known 'backward' Kol- 



mogorov equation (Gardiner, 2004) in terms of its initial condition: 



d 
dt 



A{t\xo) 



A{t\xo) 



(40) 



We can now integrate the above equation over all times to obtain an ordinary 
differential equation for T (xq). Initially, ^(0|a;o) = 1 — e~^° so that: 



^(oo|xo) — ^(0|xo) = e" 



XO 



T(xo) . 



(41) 



We note the following boundary condition on the extinction time: T (0) = 
0. This simply means that a population starting with a zero Poisson mean 
is extinct immediately, as can also be verified from the absorbing boundary 
condition on A. This condition can be utilized to solve the extinction time 



differential equation, Eq (41) 



We first introduce an auxiliary function. 



Xo {b — cxq 



dXO 



2A{x] 
D{x) 

XO 



exp 



Xo 



{b - cxoY 



dx 



1 - — — 1 dx 
b — cx 



(42) 
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where z/ = ajc— 1. In terms of?/', Eq (41) has the particular solution: 



ax 







This satisfies the backward Kolmogorov equation, and the boundary condition 
at a; = 0. In principle, another solution is possible, as the differential equation 



is of second order. However, Eq (41) has a regular singular point at x = Xm, 
whose indicial equation indicates that any other independent solution of the 
homogeneous equation would be singular at d hence inadmissable. 

This result is remarkably similar to the standard expression for a first-passage 



time of a diffusion process (Gardiner, 2004) except for the factor of (1 — e ^) in 



the integrand. This has the simple intuitive interpretation that it projects out 
the fraction of the population in a given Poisson ensemble that is still 'alive', 
ie, has not yet reached extinction. While this term is essential if an exact result 
is required, it is typically a small correction in an expression dominated by 
the exponential factors in ip. 

Extinction time results for typical parameter values are given in Table (), com- 
pared to results for direct Gillespie simulations. There is agreement within two 
standard deviations of the computational sampling error in all cases. This pro- 
vides numerical evidence for the equivalence between the discrete and contin- 
uous variable techniques for calculating extinction times. Numerical integra- 



tions of Eq (43) were checked for accuracy to at least four significant figures, 
using an adaptive routine (tolerance of .5 x 10~^) for the outer integral, and 
a fixed step integration with 4 x 10^ steps of the ^/z variable to ensure good 
accuracy at small z values. 

This set of calculations and simulations investigates the mean time to extinc- 
tion for a variety of parameter combinations. The first two sets of population 
sizes are relatively small to facilitate direct comparison between the Gillespie 
method and exact approach. The Gillespie results were obtained by simula- 
tion of 10^ realizations for each of the first twelve parameter combinations. 
Computational limitations necessitated reduction of the number of Gillespie 
simulations for the last three parameter combinations, although the numerical 
integration of the exact result is straightforward. Each simulated realization 
started with a population size drawn from a Poisson distribution with mean 
of A'c = g/c. This is also the initial condition for the exact computations. 
The extinction time results versus carrying capacity at different R values are 
also graphed in Figure 3. Nearly identical results (within samphng error) were 
found in a number of stochastic simulations, using the Stratonovich form to 
give increased numerical accuracy. 



19 



Table 1 



Time to extinction in stochastic logistic population dynamics, relative to g = 1. 



R 


A T 


Exact Te 


Gillespie gTe 


Asymptotic 


///"I'll • • 1 j_* 

# Gillespie simulations 


1.2 


5 


1.820 


1.821± 0.002 


5.031 


10^ 


1.5 


5 


4.587 


4.583±0.004 


6.580 


10^ 


2.0 


5 


10.126 


10.13 ± 0.01 


11.18 


10^ 


3.5 


5 


32.91 


32.97 ± 0.03 


32.32 


10^ 


6 


5 


83.92 


83.95 ± 0.08 


79.99 


10^ 


1.2 


10 


3.996 


3.990± 0.004 


5.534 


10^ 


1.5 


10 


12.86 


12.84 ± 0.01 


11.97 


10^ 


2.0 


10 


41.22 


41.27 ± 0.04 


36.66 


10^ 


3.5 


10 


291.3 


291.5 ± 0.3 


276.9 


10^ 


6 


10 


1430 


1431 ± 1.4 


1399 


10^ 


1.2 


20 


10.60 


10.60±0.01 


9.472 


10^ 


1.5 


20 


66.03 


66.06 ±0.06 


56.09 


10^ 


2.0 


20 


593.9 


589.7±1.9 


557.6 


10^ 


3.5 


20 


2.835 X 10^ 


2.825 X 10^ ± 89 


2.874 X 10^ 


iqs 


6 


20 


5.884 X 10^ 


5.936 X 10^ ± 6.0 X 10=^ 


6.053 X 10^ 


10^ 



A clear feature of the results is the exponential increase of extinction times 
with total population number. In addition, the results show that there can be 
changes of just as large a magnitude when the reproductive ratio R of birth 
to death rates change. This is graphed directly in Figure 4. 



These effects will be analyzed in more detail in the remainder of this section, 
where we derive the approximate analytic result found in the last column of 
the table. 

4-2 Quasi steady-state 

In order to understand these results as a diffusion process, consider the effect 
of an artificial reflecting boundary at a; = £. One can imagine this as being 
caused by the external intervention of a benevolent ecologist, wishing to pre- 
vent extinction by adding further individuals when the population becomes 
too small to be viable. 

With this new boundary there is a steady-state equihbrium, which allows us 
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10 



15 
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N 



Figure 3. Dimensionless time to extinction versus carrying capacity (population) 
N, for the same range of R values as in the table. Starting from the lower lines, 
R = 1.2, 1.5, 2, 3.5, 6. Solid lines are exact from Eq(|43l), dashed line are the steepest 



Eq (531. 



descent result from Eq (51), dotted lines are asymptotic (large N) expressions from 



to derive an effective potential for the Fokker-Planck equation. The solution 
is then: 



2Af 
' D{x) 



exp 



2A{z) 



dz 



--AT- ih 

X 



cx] 



(44) 



Apart from the normahzation factor A^, the steady-state distribution is just 
the auxiliary function used to calculate extinction times given above. This 
distribution is not normalizable in the limit of e — > 0, which indicates it is 
not a true steady-state. With this isolated population model, the steady-state 
probability of having x < e becomes infinite, for arbitrarily small values of e. 
Thus, while a local equilibrium can be reached over short times which has a 
quasi-steady-state behaviour, the only true steady-state has population zero. 

While we already have a result for the extinction time, we can use the concept 
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101.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 

R 

Figure 4. Dimensionless time to extinction versus reproductive ratio for 
N = 10, 20 (lower and upper curves respectively). Solid line is the exact result, 
dashed line is the asymptotic expression (53) 

of a quasi-steady-state to develop approximate expressions which give a better 
intuitive understanding. 



^.5 Approximate extinction times 



This quasi-steady-state distribution is typically double-peaked. Some popula- 
tions are near the deterministic steady-state, while some populations are zero. 
This may be thought of as due to a potential barrier. There is a finite prob- 
ability that a population which is locally stable will penetrate the potential 
barrier and reach the irreversible state of a; = 0, which means that an extinc- 
tion has occurred. We can calculate the extinction time approximately as an 
escape probability, from the deterministic or locally stable value through to 
extinction at a; = 0. This is equivalent to a steepest descent approximation 



(Gardiner, 2004) to the exact integral expression given in Eq (43). 



To simplify this calculation, we first change to a new variable with constant 
diffusion rate. The variable change is defined so that: 

dx I 

— = A = -y/x (6 - ex) . (45) 
ay 
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The Fokker-Planck equation for the corresponding probability distribution 
g{y) = Af{x) is: 



d_ 
di 



9{y) 



dy dy 



9{y) > 



(46) 



where Viy) = dV/dy , and V{y) is the potential for the new distribution. 
This is obtained from the quasi-steady-state distribution - since 



y(y) = -ln(A/oo) 

= —X + - In X {b — cxY~'^'^^'^ 



The potential derivative, which defines the turning points, is: 



v'{y) 



cx^ — {g + c)x + 



(47) 



(48) 



The potential has turning points at V'{y^) = , that are given by solving the 
corresponding quadratic in the Poisson variable x, so that: 



X 



2c 



(49) 



where 7 = ^J{g -\- c)^ — 26c. We now follow standard techniques to calculate 
first-passage times, and neglect the small Poisson correction of (1 — e^"^) in 
the integrand. It is also necessary to calculate the curvature of the potential, 
which is: 



=±7. 



(50) 



The average extinction time is governed by the potential difference and cur- 
vature. Provided g ^ \/2hc , it is given approximately by: 



T^^exp 
7 



V{x-) - V{x^ 



(51) 



4-4 Limiting behavior 



The table and graphs show an exponential increase in extinction times with 
carrying capacity as well as a further marked dependence on the ratio of 
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gain to death rate. Here R = b/a is the birth-death ratio, also called the re- 
productive ratio, which we introduce to obtain a dimensionless expression. We 
can make a further simplification by calculating the limit of the approximate 
exponential term at large carrying capacity, or c — >^ . This gives the limiting 
result for the turning points that: 



X 



X 



■■b/2g 



+ 1 



b 

2g 



(52) 



Next, using these asymptotic values, and evaluating the potentials at the turn- 
ing points, we find to leading order that: 



T = To exp {N, [1 - In (R) / (R - I)]} + O (1/iVj 

= roexp{iVe//} . (53) 

Here we define a fundamental time-scale of Tq, and an effective carrying ca- 
pacity Neff which determines the exponent: 



nn 27ri? 1 



ge V 2N{R~1) 

N,ff = N,[l-ln{R)/{R-l)] . (54) 

This expression is obtained assuming that R -C Nc, since we are interested 
here in the limit of large carrying capacity. It is instructive to consider what 
happens at large and small relative death rates, which gives a leading order 
asymptotic results of: 



^im ^ In T = N^jf _ liV, (i? - 1) + O (In iVj 

Hm In T = N^jf ^N^ + O (In N^) . (55) 



A comparison of the approximate method with exact results is given in Ta- 
ble ([T]). The method performs worst for the first set of population parame- 
ters, which exhibit very rapid extinction rates. This is due to the fact that 
many populations in this regime are below the critical threshold immediately, 
so don't have a quasi-steady-state. However the approximation behaves very 
well for larger population sizes where a genuine quasi-steady-state exists as 
demonstrated by the Nc = 20 simulations, where the approximation is gener- 
ally within 10% of the exact results. 
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Figure 5. Relative error in time to extinction, iTapprox — T) /T versus carrying capac- 
ity (population) A^, for R=1.5 (black line), R=2.0 (blue line) and R=3.5 (green line) 
. Small R values are the least accurate (uppermost) lines at the = 5 intercepts. 

This approximate result is also compared to the exact calculation in Figure 5, 
where the relative error is graphed for three different R values. Agreement is 
excellent for large populations, provided that growth is larger than critical, ie, 
9 ^ 9c = V2bc , which means that Nc 3> 2R/{R — 1). If the growth rate is less 
than this critical value, the potential has no local minimum, and extinction 
will occur on time-scales of the order of the inverse growth-rate for any initial 
population. We found no advantage in using the full steepest descent result 



of Eq (51) over the simpler equation (53), with the asymptotic result actually 



giving smaller errors at low values. 



In summary, the asymptotic results from the steepest-descent method are 
accurate to within a few percent at large carrying capacity, but are not reliable 
below the critical carrying capacity, where exact results or direct simulations 
can be used. 

These asymptotic results dramatically show that while large carrying capacity 
increases extinction times exponentially, it is not the only factor involved. 
There is also an exponential dependence on the reproductive ratio. For a given 
growth rate and carrying capacity, populations are exponentially shorter-lived 
at small reproductive ratio, as i? — > 1 . 

In summary, with carrying capacity above the critical value, there is a quasi- 
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stable population at ,x'+ ~ g/c. A population of this size can last an exponen- 
tially long time, although it is not globally stable. There is a saddle-point, or 
local maximum in the potential at a minimum critical population oi x = x~ . 
If fluctuations occur below this minimum critical population of ~ b/2g, the 
population is too low to be sustainable. This means that extinction becomes 
likely over a short time-scale. Extinction is exponentially much more rapid 
with small carrying capacity Nc, and also with a small reproductive ratio i?, 
which leads to large demographic fluctuations. 



5 Stochastic discrete logistic model 

So far we have only included demographic noise, due to the intrinsic random 
behavior of discrete jump events. In addition to this, there are external fluc- 
tuations, due to variations in environmental parameters like temperature that 
affect food supply or other factors relevant to survival. In our approach, these 
are represented as random, ie stochastic, time-dependent rate constants k{t). 
In general, knowledge of these rate constants and their statistical fluctuations 
would allow more accurate predictions of average extinction times. We know 
that the rates have a finite correlation time; no changes are instantaneous. 
Despite this, it is useful to treat the limit of a short correlation time relative 
to the average growth rate g. It is this limit that we consider here due to its 
analytic simplicity. More general results are possible that include the effects 
of finite correlation times, but wiU not be included here. 

Effects of this type are sometimes treated by including fluctuating terms at 
the rate equation level, without demographic noise due to discrete events. This 
has the drawback that the logistic model with time-varying rates cannot lead 
to extinction. It simply has no absorbing state. 

This problem may be circumvented approximately by assuming that a small 
population - say a; = 1 - is equivalent to extinction. However, as we have 
shown in the previous section, demographic noise itself plays a role in causing 
fluctuations leading to extinctions, even in a static environment. Further, the 
critical population leading to a high likelihood of demographic extinction is 
not necessarily at x = \, but instead is at x^ ~ b/2g. This depends critically 
on the ratio of birth to growth rates. Assuming that a; = 1 is equivalent to 
extinction can therefore lead to a serious over-estimation of extinction times 
if death and birth rates nearly cancel. 

In the following section we develop a unified theory that includes both envi- 
ronmental and demographic fluctuation effects in calculating the extinction 
rate. 
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5.1 Environmental noise 



We now include a specific model for the environmental noise. Environmen- 
tal parameters are the time-dependant rates in the logistic equation, which 
are positive, and physically have finite bandwidth fluctuations. These can be 
modeled by writing the environmental rate parameters as ki{t) = exp(0j(t)). 
Here (pi describes the fluctuations, and is modeled using a stochastic differen- 
tial equation of form: 

dcj) = [(t)dt + (TiciW^^™] , (56) 



where dWi is an external environmental noise-source such that {dWidWj) = 
Sijdt, and 7^, (Tj describe the rate of change and the relative noise variance of 
the i — th rate parameter. The resulting master equation with time-dependent 
rates can be termed the stochastic discrete logistic equation, or SDLE. 

For analysis, the time-varying rates can simply be inserted directly into the 



demographic Stratonovich equation, (31). This can be numerically integrated 



with any model of environmental stochasticity. 

If the environmental fluctuation time-scales 7"^ are much smaller than the de- 
mographic time-scales g~^, adiabatic ehmination wfll result in an approximate 
broad-band stochastic noise with variance {k^aiY5{t). Provided this variance 
is relatively small compared to the rate, the fluctuations can then be linearized 
to give additive environmental noise in the rates. 

Whfle one can solve the stochastic equations numerically without these sim- 
plifications, the result is more readily treated analytically, and wfll be treated 
in detail in the remainder of this section. As a simple example of this, we 
wfll consider a fluctuating death rate with environmental noise variance in the 
broad-band limit given by: 

{5a{t)5a{t')) ~ 2a5{t - t') . (57) 



We emphasize that broad-band external noise of this type must be included in 
the Stratonovich form of the demographic stochastic equation. Unhke the Ito 
form, the Stratonovich equation has a well-deflned broad-band limit with an 
external noise source. The relative effect of the environmental noise depends on 
the competition term c, so we wifl define r = cr/c as the relative environmental 
noise. 
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5.2 Unified Stratonovich equations 



Adding a fluctuating death rate to the demographic equation in Stratonovich 



form, Eq (31), we obtain: 



[g + c- 6a{t)] X - b/2 



cx 



dt + J2 {hx - cx^)dW 



(58) 



Because the demographic and environmental fluctuations are assumed to be 
independent, their variances can be added. As both competition and environ- 
mental noise are now broad-band noise sources proportional to x, we introduce 
c to describe the combined noise variance, where 

c = c — a = c(l — r) . (59) 



It would not be correct to make this change in (29), the demographic Ito 
stochastic equation; this does not correspond to the broad-band limit of a 
finite bandwidth stochastic process. 

Environmental noise has the opposite effect to demographic noise from density- 
dependent competition. It increases the variance to super-Poissonian levels, 
whfle intra-specific competition reduces the variance. The resulting unified 
Stratonovich equation with demographic and environmental noise is: 



dx^ = ([g + c]x- b/2 - cx^) dt + ^2 {bx - dx^)dW , (60) 



As previously, we assume absorbing boundary conditions at a; = . This means 
that a negative value of x will not occur. For positive c, there is a refiecting 
boundary at x = b/c. Since x < b/c, the equation always remains real. If c is 
negative, the upper boundary is at infinity. 

A typical simulation of these equations is given in Figure 6. The dimensionless 
parameter values are r = 0.75, = 10, R = 2. These simulations use an RK4 
(Runge-Kutta) algorithm, with step-size At = 0.08. There are 1000 trajecto- 
ries in the ensemble averages. The mean extinction time, using the method 



of Eq (39), was calculated to be T = 29.18 ± 0.9 , with time-scales chosen so 
that g = I. 



5.3 Extinction times 



We now wish to analyse the combined effects of demographic and environ- 
mental noise on extinction times. This is obtained by transforming the unified 
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t 

Figure 6. Graph of Poisson population versus time: mean values given by the dotted 
line, sample simulation given by the solid line, including environmental noise with 
r = 0.75, A^ = 10, i? = 2. 



Stratonovich equation (60) back into an Ito stochastic equation using Eq (12) 



and hence to a Fokker-Planck equation of form: 



d_ 
di 



fix) 



d_ 
dx 
d_ 
dx 



cx — g 



- + — 
cx) dx 



-Aix) + ^D{x) 



x{h — cx) f (x) 
fix) . 



(61) 



We note from Eq (12), that the transformation introduces an additional term 



in the drift that originates from the environmental noise term in the Stratonovich 
equation. The quasi-stationary (steady-state) solution is then: 



ip (x) 



X {h 



cx] 



exp 



LO 



h — cx 



dx 



(62) 



Here we have introduced a new parameter a = b/{l — r) — cr — g. This equation 
can be integrated to give: 



^ (x) = (6 - Ixf^"-^ . (63) 



X 
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As with pure demographic fluctuations, an exact extinction time result is 
obtainable using: 



T{xo) = 2 



I 



dx 



j ^(^) (l - e-') dz 



(64) 







D{x)ip{x) 



X 



We note that for a > c, one formally has x^ = oo in this treatment. This 
is an artifact of the use of a Gaussian noise source, in which there is a small 
probability of an anomalous, negative death rate. However, as the large N tails 
of the distribution have little or no effect on the extinction time, we expect 
this to be relatively unimportant. It should be remarked that the detailed 
behaviour of the tails of distribution of death rate for large positive values 
are much more signiflcant, especially if there are appreciable departures from 
Gaussian statistics. Applying this result to the parameters given in Figure 
6, wc obtain T = 28.2109, within one standard deviation of the stochastic 
simulation results. 

This excellent agreement between analytic results and simulations demon- 
strates that one does not have to use this exact theory. In fact it may be better 
not to in some cases. If the environmental noise has known statistical prop- 
erties different to those assumed here, it is preferable to use the simulations 
described above, with appropriate noise sources having the true environmental 
statistics. 



5.4 Steepest Descent Approximation 

The quasi-stcady-statc distribution is still doublc-pcakcd but with appreciably 
smaller peak potentials when a > c. Approximate results are obtained as 
before by calculating the extinction time as an escape probability, from the 
deterministic or locally stable value through to extinction at x — 0. We change 
to a new variable y, with constant diffusion rate. This variable change is now 
defined so that: 



The Fokker-Planck equation for the corresponding probability distribution 




(65) 



g{y) = Af{x) is: 



d_ 
di 



9{y) 



d_ 
dy 



(66) 
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where Vl{y) = dVa/dy , and Va{y) is the potential for the new distribution. 
This is obtained from the quasi-steady-state distribution - since 



y,(y) = -ln (A/oo) 

= —x/ (1 — r) - 



In 



x{h — cx 



l-2a/c 



(67) 



Alternatively, we can make a variable change to the Stratonovich form of the 
drift, which gives the derivative in an equivalent but simpler form as: 



V' 



1 

A 



cx"^ — {g + c)x H — 



(68) 



The potential has turning points at V^{y^) — , that are exactly the same as 
in the pure demographic case. They correspond to the points where the drift 
vanishes in the Stratonovich form of the equation, which is not changed by 
environmental noise, so that: 



X 



± 



2c 



(69) 



where 7 = y (gf — c)^ — 26c as in the purely demographic case. 

To obtain the extinction time, it is also necessary to calculate the curvature of 
the potential, which is, surprisingly, unchanged from the purely demographic 
case: 



K."(y±)=A^K:(y±) 
= ±7. 

The average extinction time is now given by: 

27r r 

7; = — exp V„{x-)-Va{x^ 

7 '- 



(70) 



This result includes both environmental and demographic contributions to the 
extinction time in a single expression. Just as in the case of pure demographic 
noise, we can simplify this expression by taking the large N hmit, while fixing 
the relative level of external noise. To simplify the resulting expressions, it is 
convenient to define an effective reproductive ratio, defined as R = b/{b+{r — 
l)g)- For a stationary environment, R = R. This, of course, is not the actual 
reproductive ratio, but rather a scaled parameter that reflects the relative size 
of fluctuations in birth events. 
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This gives rise to the following simple result, which clearly reduces to the 
purely demographic result at r = 1: 



T^Toexp 



1 — r L 
Toexp{iVe//} . 



In (^) / 



+ 0{1/N,) 



(72) 



As before, we can define a fundamental time-scale of Tq, and an effective 
carrying capacity iVe// which determines the exponent: 



Tn 



27r^V(i-r) 



1 — r 



1 — r 



g \\ 2Ne{R - 1) 
In (^) / 



1 - 



- 1 



(73) 



Just as in the previous section, one can consider what happens at large and 
small R, which gives a leading order asymptotic result of: 



Jim InT^AT,^^^— ^ +0(lniV,) 
lim In T = Neff -^Nc/{l-r) + (In N^) . 

R>1 



(74) 



Graph showing typical results are given in Figures 7 and 8 for relative environ- 
mental fluctuations of r = 1.5 and r = 3 respectively. We see that extinction 
rates for carrying capacity of = 20 are now three orders of magnitude more 
rapid than with demographic effects alone; this relative discrepancy is even 
larger at higher carrying capacity. 

A clear feature of the results is that environmental fluctuations have the great- 
est relative effect on species with a large reproductive ratio, which otherwise 
would have an extremely long extinction time in a static environment. These 
long extinction times are reduced by many orders of magnitude, even with 
relatively small environmental noise. This is because environmental noise ef- 
fects can easily be much larger than the low level of population fluctuations 
purely due to demographic causes with large R values. Since environmental 
fluctuations are practically unavoidable, we see that the exceptionally long 
lifetimes found with large reproductive ratios are probably not achievable in 
real world environments. 

The effects of varying environmental noise at fixed carrying capacity A^ are 
shown in Figure 9. This shows the strong effects of environmental noise at 
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Figure 7. Graph of scaled time to extinction versus carrying capacity (population) 
N, for the same range of R values as in the table, but including environmental noise 
with r = 1.5. Starting from the lower lines, R = 1.2,1.5,2,3.5,6. Solid lines are 



the exact results. Asymptotic results using steepest descent from Eq ( 73 1 are shown 
with dashed lines. 

large reproductive ratio. It also demonstrates that the much smaller extinction 
times found with low reproductive ratios are not as sensitive to these external 
effects. This is because small R values mean high death and birth rates, so 
that the fundamental demographic noise is high. In this situation, extinction 
is always rapid. Hence the faster extinctions due to environmental noise are 
not so dramatic, although still very significant when the noise is increased 
further. 



Discussion 



We have described a simple class of models that can be used to describe an iso- 
lated self-regulating population, including both demographic and environmen- 
tal fiuctuations. We believe the stochastic discrete logistic model represents 
an appropriate basis for a new synthesis of population genetics, mathematical 
epidemiology and theoretical ecology. We have used three equivalent tech- 
niques for analysing population dynamics and extinction times: direct master 
equation simulations, stochastic equations and exactly soluble Fokker-Planck 
equations. 



33 






10 5 



10 



15 



20 



N 



Figure 8. Graph of scaled time to extinction versus carrying capacity (population) 
N, for the same range of R values as in the table, but including environmental noise 
with r = 3. Starting from the lower lines, R = 1.2,1.5,2,3.5,6. Solid lines are the 



exact results. Approximate results using steepest descent from Eq (73) are shown 
with dashed lines. Parameters are the same as in Fig 3. 

We emphasize that in the Poisson representation used here, all three tech- 
niques are exact and give identical results. This is to be contrasted with previ- 
ous work using truncated forms of the Fokker-Planck equation, in which there 
can be exponentially large errors introduced by the diffusion approximation. 

The Fokker-Planck method has the useful feature that it allows precise yet 
analytically tractable calculations of the asymptotic extinction times for large 
carrying capacity. At large we find a general dependence of log T oc N^, but 
with a very different constants depending on the reproductive ratio R and the 
relative environmental noise r. A single analytic expression agrees to within a 
few percent of the exact extinction times over a wide range of parameter values, 
if the carrying capacity is not too low. Provided these parameter values can be 
estimated, this provides a useful basis for risk analysis of survival probabilities 
in small, isolated populations. 

Our results show that extinction times have an exponential dependence on 
carrying capacity above a critical carrying capacity that depends on the re- 
productive ratio. Surprisingly, the effect of small reproductive ratio is just as 
serious as small carrying capacity, in that both cause exponential reductions 
in extinction time to the point of almost total non- viability in the short term. 
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Figure 9. Graph of scaled time to extinction versus relative environmental noise r 
for the same range of R values as in the table, at a carrying capacity of = 20. 
Starting from the lower lines, R = 1.2, 1.5,2,3.5,6. Solid lines are the exact results. 



Approximate results using steepest descent from Eq (73) are shown with dashed 
lines. Parameters are the same as in Fig 3. 

Environmental fluctuations have a similar effect, although these are much 
more pronounced when extinction rates are low, as occurs with large R values. 
These results show clearly that carrying capacity or observed population sizes 
in isolated ecosystems are not by themselves a reliable guide to long-term 
viability. Instead, the total picture of reproductive ratio and fluctuations in 
growth rates due to environmental causes needs to be included as well in any 
risk assessment. 



There are many ways in which this work can be extended. The most obvious 
is to extend this analysis to include genetics so that extinction would rep- 
resent the loss of an aUele from the population. By introducing mutation, a 
revised neutral theory of evolution could be developed that accommodated 
demographic fluctuations directly. Simflarly natural selection could be mod- 
eled either through differential growth rates, differential death rates or inter- 
genotype competition. 

A more challenging direction of enquiry will be to develop a framework for 
inference under the spatial logistic model analogous to that which Kingman's 
coalescent provides for analyzing the idealized Wright-Fisher and Moran pop- 
ulation models. 
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Figure Captions 



Figure 1 Direct simulation results using Gillespie algorithm with param- 
eters 6 = 2, a = 1, c = 0.2 . This corresponds to A'c = 5, i? = 2, (? = 1. 
The initial condition of A'o = 8 was sampled from a Poissonian with mean of 
xo = 5. The sohd hne is a single stochastic realization, showing integer jump 
behaviour. The dotted hne is an average of 10000 realizations, showing the 
exponential decline in average population leading to extinction. 



Figure 2 Poisson simulation results using Stratonovich equations with pa- 
rameters 6 = 2, a = 1, c = 0.2 or i? = 2, A'c = 5 as in Figure 1. The initial 
condition was a Poissonian with mean of xq = 5. The solid line is a single 
stochastic reahzation, showing continuous stochastic behaviour. The dotted 
line is an average of 10000 realizations, showing the exponential decline in 
average population leading to extinction, just as in the GiUespie result. Inte- 
gration step-size utihzed was dt = 0.04, RK4 Runge-Kutta algorithm. Average 
extinction time T = 10.04 ± 0.27 using EqQ 



Figure 3 Dimensionless time to extinction versus carrying capacity (pop- 
ulation) A^, for the same range of R values as in the table. Starting from the 
lower lines, R = 1.2,1.5,2,3.5,6. Sohd hues are exact from Eq(43), dashed 
hne are the steepest descent result from Eq (51), dotted hues are asymptotic 
(large N) expressions from Eq (53). 



Figure 4 Dimensionless time to extinction versus reproductive ratio R, for 
N = 10, 20 (lower and upper curves respectively). Solid hne is the exact result, 
dashed line is the asymptotic expression ( |53) ) 



Figure 5 Relative error in time to extinction, (Tapprox—T) /T versus carrying 
capacity (population) A^, for R=1.5 (black hne), R=2.0 (blue line) and R=3.5 
(green line) . Small R values are the least accurate (uppermost) hues at the 
N = 5 intercepts. 
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Figure 6 Graph of Poisson population versus time: mean values given by the 
dotted line, sample simulation given by the solid line, including environmental 
noise with r = 0.75, N = 10, R = 2. 



Figure 7 Graph of scaled time to extinction versus carrying capacity (popu- 
lation) A^, for the same range of i? values as in the table, but including environ- 
mental noise with r = 1.5. Starting from the lower lines, R = 1.2, 1.5, 2, 3.5, 6. 
Solid lines are the exact results. Asymptotic results using steepest descent 
from Eq (73) are shown with dashed hues. 



Figure 8 Graph of scaled time to extinction versus carrying capacity (pop- 
ulation) N, for the same range of R values as in the table, but including envi- 
ronmental noise with r = 3. Starting from the lower lines, R = 1.2, 1.5, 2, 3.5, 6. 
Solid lines are the exact results. Approximate results using steepest descent 
from Eq (73) are shown with dashed hues. Parameters are the same as in Fig 
3. 



Figure 9 Graph of scaled time to extinction versus relative environmental 
noise r for the same range of R values as in the table, at a carrying capacity 
of A^ = 20. Starting from the lower lines, R = 1.2, 1.5, 2, 3.5, 6. Solid hues are 
the exact results. Approximate results using steepest descent from Eq (73) are 
shown with dashed lines. Parameters are the same as in Fig 3. 
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